#include "AnaDef.h"
#include "XGLUtils.h"
#include "AliTRDdEdxUtils.h"

void fit(TH1D *hh, TList *ll)
{
  Double_t par[3]={hh->GetEntries(), hh->GetMean(), hh->GetRMS()};
  XGLUtils::ChisquareFit(hh, XGLUtils::AGaus, 3, par);

  TF1 *ff=new TF1(Form("func%s",hh->GetName()), XGLUtils::AGaus, -2, 2, 3);
  ff->SetLineColor(kRed);
  ff->SetParameters(par);
  ff->Draw("same");
  ll->Add(hh);
  ll->Add(ff);
}

void getResolution(const Int_t itype)
{
  //========================================================
  //========================================================
  gErrorIgnoreLevel = 1001;
  TTree *tree=(TTree*)gDirectory->Get("newdraw_Tree");
  
  TList *ll=new TList;  
  ll->SetOwner();

  TCanvas *cc=0x0;
  TH1D *hh=0x0;
  TH2D *h2=0x0;
  //-
  TString basecut("ifill==0 &&");
  TString trdcut;

  //elmass: 5.10999e-04
  //pimass: 1.39570e-01  

  tree->SetAlias("rtpc","tpcsig/AliTRDdEdxUtils::QMeanTPC(pow(10,ltbg))");

  TString tpcbg("AliTRDdEdxUtils::QMeanTPC(TMath::Power(10,ltbg))");
  TString tpcmom("AliTRDdEdxUtils::QMeanTPC(mom/1.39570e-01)");

  TString trdQ0bg, trdQ0mom;
  if(itype==2){
    trdQ0bg = "AliTRDdEdxUtils::Q0MeanTRDpp(TMath::Power(10,ltbg))";
    trdQ0mom = "AliTRDdEdxUtils::Q0MeanTRDpp(mom/1.39570e-01)";

    basecut = Form("itype== %d && ", itype);
    basecut += "run >= 126461  && run <= 130930 &&"; printf("selecting 10eeeeeeeeee !!\n");
    trdcut="&& trdnch0>=6 && rtpc>0.8 && rtpc<1.2 ";
  }
  else if(itype==4){
    trdQ0bg = "AliTRDdEdxUtils::Q0MeanTRDPbPb(TMath::Power(10,ltbg))";
    trdQ0mom = "AliTRDdEdxUtils::Q0MeanTRDPbPb(mom/1.39570e-01)";

    basecut = "itype >=10 && itype <=11 && ";
    trdcut="&& trdnch0>=6 && rtpc>0.95 && rtpc<1.3 ";
  }
  else{
    printf("error itype! %d\n", itype); exit(1);
  }

  basecut += "(pid==0 || pid==2 || pid==4)";
  
  printf("basecut: %s\n", basecut.Data());
  printf("trdcut: %s\n", trdcut.Data());

  //========================================================================
  //========================================================================
  cc=new TCanvas("crtpc") ;gPad->SetGrid();gPad->SetLogz();  ll->Add(cc);
  hh=new TH1D("hrtpc","",100, 0.5, 1.5);
  tree->Draw(Form("rtpc>>%s", hh->GetName()), basecut);
  fit(hh, ll);

  cc=new TCanvas("crEltpc") ;gPad->SetGrid();gPad->SetLogz();  ll->Add(cc);
  h2=new TH2D("hrELtpc","",100,0.3, 5,100, 0, 5);
  tree->SetAlias("rELtpc",Form("tpcsig/%s", tpcmom.Data()));
  tree->Draw(Form("rELtpc:mom>>%s", h2->GetName()),basecut,"colz");
  ll->Add(h2);

  cc=new TCanvas("cdEltpc") ;gPad->SetGrid();gPad->SetLogz();  ll->Add(cc);
  h2=new TH2D("hdELtpc","",100,0.3, 5,100,-50, 50);
  tree->SetAlias("dELtpc",Form("tpcsig-%s",tpcmom.Data()));
  tree->Draw(Form("dELtpc:mom>>%s", h2->GetName()),basecut,"colz");
  ll->Add(h2);
  //------------------------------------------------------
  //------------------------------------------------------

  tree->SetAlias("trdQ0sig","trdQ0");

  cc=new TCanvas("crtrdQ0") ;gPad->SetGrid();gPad->SetLogz();  ll->Add(cc);
  hh=new TH1D("hrtrdQ0","",100, 0.5, 1.5);
  tree->SetAlias("rtrdQ0",Form("trdQ0sig/%s", trdQ0bg.Data()));
  tree->Draw(Form("rtrdQ0>>%s", hh->GetName()),basecut+trdcut);
  fit(hh, ll);

  //------------------------------------------------------
  //------------------------------------------------------
  const Int_t nx = 100;
  const Double_t xmin = 0.3;
  const Double_t xmax = 5;
  const Int_t ny = 100;
  const Double_t ymin = 0;
  const Double_t ymax = 2;

  cc=new TCanvas("crEltrdQ0") ;gPad->SetGrid();gPad->SetLogz();  ll->Add(cc);
  h2=new TH2D("hrELtrdQ0","",nx, xmin, xmax,ny, ymin, ymax);
  //tree->SetAlias("rELtrdQ0",Form("trdQ0sig/%s", trdQ0mom.Data()));
  //tree->Draw(Form("rELtrdQ0:mom>>%s", h2->GetName()),basecut+trdcut,"colz");
  tree->Draw(Form("trdQ0sig:mom>>%s", h2->GetName()),basecut+trdcut,"colz");
  ll->Add(h2);

  const Double_t pimass = 1.395700e-01;
  const Double_t elmass = 5.109989e-04;
  const Double_t prmass = 9.382720e-01;
  TH1D *mpi=new TH1D("mpi","",nx, xmin, xmax);
  TH1D *mel=new TH1D("mel","",nx, xmin, xmax);
  TH1D *mpr=new TH1D("mpr","",nx, xmin, xmax);
  for(Int_t ii=1; ii<=mpi->GetNbinsX(); ii++){
    const Double_t xx= mpi->GetXaxis()->GetBinCenter(ii);
    const Double_t piy = itype==2 ? AliTRDdEdxUtils::Q0MeanTRDpp(xx/pimass): AliTRDdEdxUtils::Q0MeanTRDPbPb(xx/pimass);
    const Double_t ely = itype==2 ? AliTRDdEdxUtils::Q0MeanTRDpp(xx/elmass): AliTRDdEdxUtils::Q0MeanTRDPbPb(xx/elmass);
    const Double_t pry = itype==2 ? AliTRDdEdxUtils::Q0MeanTRDpp(xx/prmass): AliTRDdEdxUtils::Q0MeanTRDPbPb(xx/prmass);
    mpi->SetBinContent(ii,piy);
    mel->SetBinContent(ii,ely);
    mpr->SetBinContent(ii,pry);
  }
  TH1D *hhs[]={mpi, mel, mpr};
  for(Int_t ii=0; ii<3; ii++){
    hhs[ii]->SetLineWidth(2);
    hhs[ii]->Draw("same C");
    ll->Add(hhs[ii]);
  }

  //======================

  TFile *fout=new TFile(Form("resolution%d.root",itype),"recreate");
  ll->Write();
  fout->Save();
  fout->Close();
}

int main(int argc, char *argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=3){
    printf("argc!=3\n");
    return 1;
  }

  TFile::Open(argv[1]);

  getResolution(atoi(argv[2]));
  return 0;
}
